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ABSTRACT 

I describe a nested-grid particle-mesh (NGPM) code designed to study gravita- 
tional instability in three-dimensions. The code is based upon a standard PM code. 
Within the parent grid I am able to define smaller sub-grids allowing us to substan- 
tially extend the dynamical range in mass and length. I treat the fields on the parent 
grid as background fields and utilize a one-way interactive meshing. Waves on the 
coarse parent grid are allowed to enter and exit the subgrid, but waves from the sub- 
grid are precluded from effecting the dynamics of the parent grid. On the parent grid 
the potential is computed using a standard multiple Fourier transform technique. On 
the subgrid I use a Fourier transform technique to compute the subgrid potential at 
high resolution. I impose quasi-isolated boundary conditions on the subgrid using the 
standard method for generating isolated boundary conditions, but rather than using 
the isolated Green function I use the Ewald method to compute a Green function on 
the subgrid which possesses the full periodicity of the parent grid. I present a detailed 
discussion of my methodology and a series of code tests. 

Key words: cosmology - galaxies: clustering - numerical methods. 



1 INTRODUCTION 

Over the past decade N-body techniques have become the 
dominant method for studying the clustering of mass on 
large scales (see Efstathiou et al. 1985, or Hockney & East- 
wood 1981). Direct particle-particle (PP) codes were the 
first n-body codes which were widely used, and they reached 
a remarkable state of development as discussed by Aarseth 
(1985). These codes tend to be very time consuming since 
each particle interacts with every other particle directly via 
a 1/r 2 force. This technique allows for very high spatial res- 
olution but at the expense of CPU time. As the particles 
clump the the small scale dynamics dominate and the CPU 
time required jumps dramatically. Treecodes are dramati- 
cally faster since they compute direct PP interactions only 
locally and the far field computations are performed using 
nodes of particles. These codes do not require a mesh and 
hence do not have the spatial resolution problems associated 
with the introduction of a grid. For both of these codes the 
small-scale spatial resolution is set by the finite softening 
length which is introduced to avoid the formation of tight 
binary pairs which would force a significant decrease in the 
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time step and lead to dramatic CPU requirements. Particle- 
Mesh codes have circumvented a number of those problems 
at the expense of the finite spatial resolution introduced by 
the existence of the grid system. The forces are interpolated 
from the grid and hence they are truncated on small scales. 
This helps to eliminate the two-body effects present in PP 
codes. Thus PM codes are well suited to the study of Vlasov- 
like systems where it is essential to suppress two-body ef- 
fects, such as the dark matter density field in cosmological 
simulations. In an attempt to increase the spatial resolution 
of PM codes, Particle-Particle-Particle-Mesh (P 3 M) codes 
were developed. P 3 M codes are a hybrid class which use PM 
methods to compute large-scale forces and PP methods for 
small-scale interactions. This modification to PM codes par- 
tially helps to increase resolution. On the other hand they 
will suffer the same dramatic slow down which occurs with 
PP codes as clumping becomes significant. Furthermore, as 
emphasized by Sellwood (1987) P 3 M codes may introduce 
two-body effects on small-scales thus making them unsuit- 
able for use in modeling Vlasov-like systems. Pen (1994) has 
constructed a code which he refers to as a "Linear Moving 
Adaptive PM Code" . This code is adaptive in the traditional 
sense, where the mesh spacing varies according to some local 
quantity, in his case the density. As a consequence of Adap- 
tive Mesh Refinement (AMR) Pen finds that the particles 
in his code feel a self-force, this could be a problem for some 
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types of initial conditions. Katz & White (1993) introduced 
a "multi-mass" technique which is based upon the hybrid 
N-body/hydrodynamics code TREESPH, which was devel- 
oped by Hernquist & Katz (1989). Katz & White (1993) 
used their "multi-mass" code to examine the properties of 
simulated galaxy clusters. In this method they use a series 
of nested lattices, with the particle mass growing smaller as 
one moves to more finely spaced lattices. Thus, they are able 
to obtain simultaneously high force accuracy and high mass 
resolution. This method is similar to the nested-grid meth- 
ods that I will discuss next. In a similar spirit to P 3 M codes 
where one increases the force resolution without concern for 
increasing mass resolution Suisalu & Saar (1995) and inde- 
pendently Jessop, et al. (1994) have developed multi-grid 
based adaptive particle codes. Both codes are capable of 
adaptively modifying the underlying grid to accommodate 
the varying distribution of particles. Suisalu & Saar (1995) 
use a full multi-grid method for the potential solver which 
is able to refine on regions of high density increasing their 
force resolution. Jessop, et al. (1994) use a relaxation Poisson 
solver and Neumann boundary conditions interpolated from 
the parent grid, coupled with a methodology for adaptively 
creating sub-grids in regions where increased force resolution 
would be of benefit. 

Nested mesh codes have been available in the atmo- 
spheric sciences for a number of years, and Koch & Mc- 
Queen (1987) provide a nice introduction. In these fields it 
was realized that nesting a fine mesh within a coarse mesh 
can be a very economical way to achieve higher resolution 
without increasing dramatically the required memory and 
CPU resources that an equivalent larger grid would need. 
Nested mesh techniques are one way to extend resolution 
in simulations of large-scale structure. Peebles (1980) first 
argued that long wavelength Fourier modes can couple to 
much shorter wavelength modes increasing the power on 
small scales. More recently Jain & Bertschinger (1993) have 
come to the same conclusion. As a consequence attempts to 
increase resolution by reducing the box size appear to be 
doomed from the start. Furthermore, moving from simple 
PM codes to P 3 M codes on a similar sized grid appears to 
be capable of increasing the resolution by only a small fac- 
tor, roughly a factor O(10 3 ) shy of the required resolution 
for full simulations capable of addressing both the details of 
galaxy formation while still being able to follow the large 
wavelength modes. 

Over the last several years a number of nested grid 
codes have been developed for use in cosmology and/or as- 
trophysics. Chan et al. (1986) developed a nested mesh code 
for use in studying galaxy collisions. In the code of Chan, 
et al. the sub-grid potential is computed using a standard 
finite difference approach and the boundary conditions for 
the sub-grid are interpolated from the known coarse grid po- 
tential. Villumsen (1987) constructed a Hierarchical Particle 
Mesh (HPM) code which is similar to the code I will describe 
here. Recently, Anninos, et al. (1993) have reported a nested 
mesh code for use in cosmology which not only evolves the 
collisionless dark matter on a fine grid but performs a hydro- 
dynamic calculation on the fine grid as well allowing them to 
follow both the dark matter and the baryonic matter compo- 
nent. The code developed by Anninos et al. is similar to that 
developed by Chan et al. in the sense that boundary condi- 
tions for the sub-grid region are obtained directly from the 



parent grid by interpolating potential values on parent mesh 
cells to the sub-grid. In a related approach to the nested grid 
schemes discussed above Couchman (1991) has modified the 
standard (P 3 M) algorithm by introducing refined meshes in 
regions of high density. The traditional complaint against 
such codes has been that as clustering evolves they tend to 
slow down because there are an increasing number of par- 
ticles within one cell of one another. This causes the PP 
step to become the dominant portion of the code and sig- 
nificantly slows down the code as a whole. To circumvent 
this, in regions of high density Couchman introduces refined 
meshes to guarantee that there is never a very large number 
of particles within the same grid cell. This prevents the PP 
portion of the code from dominating the overall runtime. 

This implementation of a nested-grid algorithm has a 
number of advantages over the methods previously men- 
tioned. Unlike codes such as P 3 M, and tree codes I am able 
to get very high force resolution but without the slowdown 
associated in those codes for highly clustered distributions. 
Furthermore, I am able to not only increase the force reso- 
lution, but also the mass resolution on the sub-grid region. 
This is something which only other nested-grid codes are at 
present capable of doing. In comparison to the nested-grid 
codes, I feel that my mass advection scheme should make 
my method easier to generalize to a more adaptive scheme. 
I also enforce a Courant-Friedrichs-Levy (CFL) condition on 
the sub-grid time step, guaranteeing that the time integra- 
tor does not go unstable. This is contrary to the Villumsen 
(1988) version of a nested-grid code where no CFL condition 
is enforced on the sub-grid particles. 

This paper is organized in the following way: section §11 
discusses in detail the algorithm that I have implemented, 
Section §111 presents the results of a number of tests of the 
code, and finally in section §IV I present my conclusions. 



2 GENERAL ASPECTS OF THE CODE 

The general philosophy of the code is quite simple. I am free 
to define a sub-grid anywhere within the parent grid. As 
parent particles move into the sub-grid region I view them 
as consisting of n 3 9 sub-grid particles which initially are on 
a cubic lattice the size of one parent grid cell. As the par- 
ent particle enters the sub-grid the sub-grid particles are 
free to move on a sub-grid mesh which is typically much 
finer than the parent mesh. This allows us to gain signif- 
icantly better length resolution, albeit in a small localized 
region. Furthermore, since the typical sub-grid particle mass 
is m S g = rripg/nig, where m pg is the parent particle mass 
I am able to gain additional mass resolution beyond that 
possible with the original PM code. 

Up to this point the discussion has focused primarily 
on methods to increase force resolution. As emphasized by 
Melott & Shandarin 1989, and Melott 1990, increasing the 
force resolution is not enough in cosmological simulations. 
The authors argue that increasing the force resolution with- 
out a corresponding increase in the mass resolution will lead 
to the growth of spurious perturbations due to shot noise 
from the particles. As a consequence any method for increas- 
ing the force resolution without increasing the mass resolu- 
tion will tend to lose phase information, and thus is not 
getting the physics correct in the model. Unfortunately, this 



© 0000 RAS, MNRAS 000, 000-000 



Nested Grid Code for Cosmological Simulations 3 



Figure 1. An example of the adjacent mesh structure that the 
nested mesh code utilizes. 

error will not show up in the autocorrelation under most cir- 
cumstances (Melott 1990). Therefore, results may very well 
be incorrect, but as long as the autocorrelation is being used 
to study the system the errors may not be apparent. One 
great advantage of a nested-grid scheme over other methods 
which try to increase the force resolution without increasing 
their mass resolution such as tree codes or P 3 M, is that a 
nested-grid code will also increase the mass resolution since 
the sub-grid particles will always have a smaller mass than 
the parent grid particles. A nested-grid code will be capa- 
ble of increasing force and mass resolution simultaneously, 
thus avoiding the introduction of any phase errors due to 
the growth of spurious perturbations in the simulation. 

2.1 Grid Layout and Parent Grid/Sub-grid 
Interfaces 

I have adopted an adjacent mesh structure, as discussed by 
Koch &: McQueen (1987), for the sub-grid/parent grid hier- 
archy. A two-dimensional example is show in Fig. 1. I am free 
to define the size of the sub-grid in terms of parent cell sizes, 
once the coordinates of the lower left hand corner of the sub- 
grid have been stated. One can use an arbitrary number of 
parent grid cells and sub-grid cells, with a uniform spacing 
on both grids. In addition, the sub-grid region is fixed and 
is not allowed to move as the simulation progresses. 

In my code a one-way interface is utilized so that the 
parent grid particles can influence the sub-grid particles, but 
the sub-grid distribution cannot back-react with the parent 
particles. Anninos et al. (1993) found that significant noise 
was generated when a two-way interface was used. As dis- 
cussed by Koch & McQueen (1987) this may not be unrea- 
sonable to expect since the higher frequency waves present 
on the sub-grid may generate false waves at the interface, 
or since the high frequency waves cannot be represented on 
the much coarser parent grid they may lead to significant 
aliasing effects. 



As emphasized by Koch & McQueen (1987) and An- 
ninos et al. it is necessary to include a buffer zone around 
the sub-grid to smooth any density or force transients which 
may be introduced as the sub-grid particles enter the sub- 
grid region. I have found that a buffer zone the width of 
two parent cells is sufficient to minimize any transient ef- 
fects from mass movement into or out of the sub-grid, but 
the width of the buffer zone is an input quantity and it can 
be varied if need be. In addition, the buffer zone is divided 
in half creating a inner and outer buffer zone. The sub-grid 
particles are moved using the parent grid potential until they 
enter the inner buffer zone. At that point the sub-grid par- 
ticles are moved using the high resolution potential of the 
sub-grid region. The use of an inner and outer buffer zone 
helps the sub-grid particles to slowly relax to the additional 
power which the sub-grid potential possesses, and therefore 
by the time the sub-grid particles have reached the actual 
sub-grid they have relaxed and are moving under the influ- 
ence of the sub-grid potential only. This helps to eliminate 
any edge effects which might come into play should the sub- 
grid particles not be allowed to relax prior to their entering 
the sub-grid region. 

Furthermore, I have found the need to add a layer of 3 
sub-grid cells in the buffer zone. Thus, increasing the mesh 
upon which the density assignment is performed by 6 grid 
cells in each direction so that the entire sub-grid is increased 
by two parent cells in all directions plus an additional 3 sub- 
grid cells in all directions. The addition of these extra sub- 
grid cells helps to guarantee that the advection of sub-grid 
particles into the sub-grid will not introduce any unwanted 
density fluctuations at the sub-grid/parent grid boundary. 
Details on how the forces are computed on particles as they 
move into the sub-grid region are given below. 

2.2 Force Calculation 

The forces on the parent particles are obtained by a simple 
two-point differencing of the potential on the parent grid. To 
obtain the forces on the sub-grid I break the calculation up 
into two stages (for a single nested mesh). First, I compute 
the density on the parent grid and set the density over the 
sub-grid region equal to the mean density on the sub-grid. 
By setting the density equal to the mean over the sub-grid 
region I erase all of the fluctuations on the sub-grid, but the 
total mass present on the sub-grid is still present. This den- 
sity field is then Fourier transformed and convolved with the 
standard seven-point approximation to the Green function 
in k space, i.e. 



6 - 2cos(^f) - 2cos(^) " 2cos(^) 

where 9 is a dimensionless parameter defined to be 9 = 
AnGpSt 2 , and G is Newton's constant, p is the global mean 
density, 5t is the magnitude of the time step. As discussed 
in Potter (1973) 6 < 0.5 for stability. I have chosen 9 to be 
0.49 . 

When this potential is inverse Fourier transformed it 
represents the contribution to the subgrid potential from 
those parent particles which have not entered the subgrid. 
I then difference this potential and interpolate using a CIC 
scheme to compute the background forces on the subgrid 
particles. This leaves us with the final set of forces, the forces 
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obtain the sub-grid contribution to the forces on the subgrid 
particles. 

Using the derivation of Ziman (1972) or Sangster & 
Dixon (1976) the Green function can be computed using 
the Ewald method. It takes the following form, 



G(x, nL) = 

V 71 " 



1 \ -> 2-7T r* _, 

+ — - > cos — h ■ x 
nL ^ L 



rh 2 /Lo 



K 2 



+ i =: erfcc(a|:r — nL\) 

^ x — nJ\ 



(2) 



Figure 2. A plot of the relative error between the Ewald gener- 
ated Green function and the isolated Green function. 



on the subgrid particles due to themselves. Here I use a 
transform technique similar to Villumsen (1989), but rather 
than imposing isolated boundary conditions on the subgrid 
potential I use "quasi-isolated" boundary conditions. 

To test the relative error that would be made by us- 
ing purely isolated boundary conditions as opposed to the 
"quasi-isolated" boundary conditions which I use in this 
code, I compute the relative error between the Ewald gen- 
erated Green function and the isolated Green function. The 
horizontal axis is the quantity r/L where r is the radius and 
L the box size. The vertical axis is the relative difference 
between the Ewald Green function and the isolated Green 
function. Therefore for a sub-grid size which is 1/4 of the 
parent grid the the largest radius is 0.25\/3 and the error as- 
sociated with using the isolated Green function rather than 
the Ewald Green function is roughly 40%. This plot is shown 
in Fig 2. 

Eastwood & Browning (1979) and Hockney & Eastwood 
(1981) provide a method for imposing isolated boundary 
conditions on the a cubical mesh. I impose isolated bound- 
ary conditions on the sub-grid using the method discussed 
by the above authors. This consists of doubling the density 
field, padding the extra regions with zeros, and convolving 
it with a Green function which has an imposed periodicity 
of the main grid length. The Green function on the sub- 
grid is computed using the Ewald method (1921) which pos- 
sesses the full periodicity of the parent grid. In this manner I 
take into account the images of the sub-grid particles which 
would not exist if I used the Green function appropriate for 
an isolated system. This is what I mean by "quasi-isolated" 
boundary conditions. The Green function is computed once 
at the beginning of the simulation, transformed and saved. 
The transformed density is convolved with the Green func- 
tion obtained from the Ewald method and Fourier trans- 
formed back to give us the high resolution subgrid potential. 
This potential can then be differenced and interpolated to 



where k — 2nh/L, and h is an integer triplet, and a is an 
arbitrary parameter introduced by the Ewald method and 
represents the length scale at which the real and k space 
expansions become important. Hernquist et al. (1991) found 
adequate convergence when a = 2/L, \h\ 2 < 10, and \x — 
nL\ < 3.6i. I adopt their values in this code. 

2.3 Time Integration 

The time integration on both the subgrid and the parent 
grid is performed using a standard leap-frog technique. The 
advantages of using the leap-frog method are discussed in 
Hockney & Eastwood (1981) and Potter (1973). In order to 
ensure stability of the leap-frog integration I use separate 
time steps for the parent grid and the subgrid. I discuss 
this choice in more detail below. On the parent grid the 
time step is constrained first by a Courant-Friedrichs-Lewy 
(CFL) condition, 



5x 



Stcfi = 0.25 x 
then secondly by a free fall time constraint 
8tff = 



II 



(3) 



(4) 



The parent grid time step is then defined to be the minimum 
of the above two estimates St pg — min(5t c fi, Stf f). On the 
subgrid the same two quantities are computed and the min- 
imum of the two chosen, St 3g . Then I compute their ratio 



and truncate it. This is the number of subgrid steps, 

N S te P S = mt(Stpg/St ag ), 



(•') 



I am allowed to take before I must take another parent step 
to update the parent grid potential. Then the actual subgrid 
time step size is chosen to be St sg = St pg /N s te P s- 

The use of asynchronous time steps for the sub-grid par- 
ticles is an improvement over using the same time step for 
both grids. The reasons for that are straightforward. I have 
chosen to use the leap-frog integration scheme to perform 
the time integrations, therefore I must enforce a CFL con- 
dition on both the sub-grid and the parent grid. This is to 
guarantee the stability of the time integration. At this point 
several approaches could be used. The first would be to com- 
pute the most restrictive time step on both the sub-grid and 
the parent grid from the CFL conditions on each, and evolve 
both meshes using the same time step. This would lead to 
large overhead since the sub-grid time step will nearly always 
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be much smaller than the parent time step, thus forcing me 
to evolve the parent grid particles for many more time steps 
than is needed for an accurate solution. Secondly, I could 
use asynchronous time steps, updating the parent particles 
only when necessary. I can justify the use of an asynchronous 
time step by the following two arguments. First the gravita- 
tional infall timescale should be an estimate of the dynam- 
ical timescale for the system. Since the gravitational infall 
timescale is proportional to the inverse of the square root of 
the density, and since I have much higher mass resolution 
on the sub-grid I expect much smaller infall timescales for 
the sub-grid, in other words the sub-grid particles should be 
evolving on a smaller timescale. Thus to accurately follow 
their motion I need a smaller time step. A second justifica- 
tion for the use of asynchronous time steps is also possible 
using several recent approximation methods. Based upon 
the work of Brainerd et a. (1993) and Bagla & Padman- 
abhan (1994) on the 'Frozen Potential' formalism we know 
that the potential does not vary too much over the course 
of evolving the model. Furthermore, Melott et al. (1995) ar- 
gue in their 'Stepwise Smoothed Potential' formalism that 
gravitational clustering is equivalent to smoothing the ini- 
tial potential on increasing scales. As a consequence objects 
which are distant and have been smoothed will tend to have 
a small effect on the evolution of the sub-grid distribution 
of particles. Thus, I should be able to evolve the sub-grid 
distribution separately from the parent grid safely. 

2.4 Mass Advection 

Crucial to any nested mesh code is how one deals with the 
advection of mass into the subgrid region. The technique 
that I have chosen to use is significantly different than ei- 
ther Villumsen (1987) or Anninos et al. (1993) who perform 
a background run with only parent particles, tag the parent 
particles which actually enter the subgrid then during a sub- 
sequent run they evolve the distribution of parent particles 
and the subgrid particles whose parents will eventually en- 
ter the subgrid. Prior to the subgrid particles entering the 
subgrid they are moved using the background potential used 
to move the parent particles. Then upon entering the sub- 
grid the finer subgrid potential with its higher Fourier com- 
ponents is used to move the subgrid particles. A possible 
disadvantage of this method is that one must evolve the full 
distribution of sub-grid particles, both those particles on the 
sub-grid and those that have not yet entered the sub-grid. 
I believe that my scheme is significantly better, because I 
avoid the overhead entailed by evolving those sub-grid par- 
ticles not on the sub-grid is nonexistent, and since the sub- 
grid particles do not exist until they are laid down on the 
sub-grid region my method will make it significantly easier 
to generalize my code to a fully adaptive setting where the 
code defines local sub-grid regions on its own for refinement 
and subsequently lays down the sub-grid particles without 
the need for running the code multiple times refining in a 
local region each runtime. 

The method I use for laying down the sub-grid parti- 
cles is somewhat similar to how particles are smeared out in 
a anisotropic manner in adaptive or tensor smooth particle 
hydrodynamics algorithms ( Martel et al. 1994, or Benz & 
Davies 1993). As the parent particles approach the buffer 
zone I consider their volume as having been smeared out 



in space, with the cloud for a parent particle being given 
by half the distance to the nearest parent particle in that 
particular direction. As a consequence as far as the sub-grid 
particle advection algorithm is concerned the parent par- 
ticles no longer have a volume equal to a parent grid cell 
cubed. In essence I am distorting the volume of the sub-grid 
particle cloud in phase space to be consistent with the phase 
space distortions of the parent particles. Once the charac- 
teristic size of the parent particle is known the coordinates 
of the sub-grid particles can be computed simply, assuming 
that they are laid down uniformly in each direction. 

The x, y, and z components of the velocity are then as- 
signed by performing a multivariate interpolation between 
the nearest 3 3 parent particles in phase space in each di- 
rection. The interpolation algorithm I use was developed 
by Renko (1988). It is a multivariate interpolation scheme 
which is capable of performing the interpolation on scat- 
tered data points. Then the subgrid particles are moved us- 
ing only the parent potential until they reach the edge of the 
buffer zone discussed above at which point they are moved 
using the high-resolution potential generated on the sub- 
grid. This should be adequate for most problems of interest, 
and the following tests demonstrate that little significant 
error is introduced in the evolved sub-grid particle distribu- 
tion. Furthermore, I believe the strength of my technique 
is that it should generalize to a fully adaptive nested grid 
code somewhat easier than the Villumsen (1988) algorithm 
which requires they be laid down before evolution of the 
distribution begins. 



3 INITIAL CONDITIONS 

It has become standard place in cosmological simulations to 
generate initial conditions by using the Zel'dovich (1971) ap- 
proximation. The Zel'dovich approximation was first used to 
by Klypin & Shandarin (1983), and Efstathiou et al. (1985) 
to generate the initial conditions for the particles in a cos- 
mological simulation. Randomly generated initial particle 
positions generate shot-noise in the initial density fields on 
all scales. In models close to equilibrium at the beginning, 
the initial behavior is dominated by the collective response 
to this initial noise, thus masking the true modes of the sys- 
tem. Traditionally this noise has been reduced by turning to 
quiet starts to generate initial conditions. Byers and Crewel 
(1970) showed that arranging the particles uniformly at the 
outset suppressed the amplitude of spurious density fluctua- 
tions due to shot-noise. Further, provided there are at least 
as many particle planes as there are grid spaces in each di- 
rection Melott (1983) and Efstathiou et al. (1985) showed 
that a regular grid spacing in each coordinate direction will 
give a quiet start to a cosmological simulation on a Cartesian 
grid. 

The initial conditions for the parent particles are gen- 
erated in the standard way: the particles are laid down on 
a uniform mesh and then perturbed away from the uni- 
form mesh through the use of the Zel'dovich approximation 
(1970) which requires the initial power spectrum for the den- 
sity perturbations. The use of a uniform mesh ensures that 
the shot-noise fluctuations in the density will be minimized. 
The resulting density perturbations are consistent with the 
initial power spectrum with random phases. 
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Villumsen (1989) has discussed a method for generat- 
ing power on scales not resolvable by the parent grid. I have 
begun writing such a code and will present the tests of the 
method in an upcoming paper, and focus only on the evolu- 
tion code in this paper. 




4 CODE TESTS 



In this section I wish to present a number of tests of the 
code. 



4.1 Force/Potential Accuracy 

The first non-trivial test, compares the accuracy of the forces 
on the parent grid and on the sub-grid. One parent particle 
is laid down the parent grid with the 4 3 = 64 sub-grid par- 
ticles laid down at the same grid location on the sub-grid. 
The sub-grid is chosen to be five parent cells on a side with 
25 3 mesh cells for a refinement factor of five, and a sub- 
grid mesh size of 0.2. The peculiar choice for the number of 
sub-grid cells is made so that the 64 sub-grid particles lie 
in one sub-grid cell only and are not smeared over the eight 
nearest neighbors when the density assignment is performed. 
The density fields are computed on each grid, and the corre- 
sponding potential field is computed for each grid. Then the 
forces are computed at 2000 randomly chosen positions on 
each grid. The results are plotted in Fig. 3, the small dots 
are the forces from the parent grid, the open circles are the 
forces determined for the refinement given above, and the 
filled circles are the forces for a refinement factor of 15. I 
note several things concerning this result. First, both force 
fields have roughly the resolution expected; the parent grid 
forces go too soft at a radius of roughly 1, where I would 
have expected them to, and the sub-grid forces go too soft 
at roughly a radius of 0.2, again exactly where I would have 
expected them to. Secondly, the sub-grid forces appear to 
be much more isotropic than the parent grid forces, as mea- 
sured by the typical scatter of the points at a given radii. 
Finally the force resolution on the sub-grid appears to scale 
correctly as one increases the resolution of the code. 

4.2 Plane Wave Collapse - The Zel'dovich 
Approximation 

The simplest and most straightforward test is a one dimen- 
sional test which is useful since it possesses an analytic solu- 
tion. The analytic solution was completed by Zel'dovich in 
1970, and has come to be known as the Zel'dovich approxi- 
mation. The test has been performed with the perturbations 
aligned along the x axis similar to Efstathiou, et al. (1985) 
for a strict PM and P 3 M code. To make the test more chal- 
lenging I run the test with the perturbation moving diago- 
nally from lower corner to the diametrically opposite corner 
of the box. For this test the code was evolved to a S rm s ~ 1-0, 
where <5 rms is the RMS variation of the global density field. 
I chose to run the code only up to a 8 rms w 1.0 because the 
Zel'dovich approximation breaks down at that point, mak- 
ing a comparison to the approximation impossible. At the 
same time spurious velocities in the two orthogonal direc- 
tions were about v y , z ~ O(10~ 3 ), which, when multiplied by 
the total runtime of the simulation, corresponds to spurious 
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Figure 3. A plot of the sub-grid and parent grid force resulting 
from a single particle. The parent grid forces are given by the 
dots, the open circles give the sub-grid forces for a refinement 
of 5, the solid circles give the sub-grid forces for a refinement of 
15, and the solid line gives the exact analytical solution from the 
Ewald method. 

particle motion in the two orthogonal directions of less than 
l/30th of a subgrid cell. Fig. 4 shows a phase space plot; 
the open squares are the parent grid particles, and the solid 
dots are the sub-grid particles. To generate this plot I project 
the particle positions onto the plane which runs diagonally 
through the box. I do not include the analytic fit because 
it would obscure most of the sub-grid points. Fig. 5 shows 
the density field for that same simulation. To make this plot 
I simply plotted the density in the cells along the diagonal 
of the box. This spurious motion is unavoidable given my 
use of quasi-isolated boundary conditions on the sub-grid. 
Furthermore it contributes to the slight appearance of noise 
that is seen in the phase space plot. Nonetheless, the motion 
appears to be a small perturbation on the overall collapse 
of the pancake, and thus the effect appears to be negligible. 

4.3 Spherical Infall 

For any code, comparison with known analytic results is 
paramount. Furthermore, tests which possess a symmetry 
not intrinsic to the code are of particular usefulness. In the 
spherical infall test I have both such cases. The code with its 
cubic cells may tend to favor tests which may possess planar 
symmetry such as the plane wave test above. The analytic 
solution to spherical infall was worked out in detail by Fill- 
more & Goldreich 1984, Bertschinger (1985). Bertschinger 
found for small values of the radius the solution approached 
a power law form. In this test a uniform distribution of par- 
ticles was laid down on both the parent grid and the sub- 
grid. Then I took several time steps in which the density 
field was given by the CIC values from the particles plus 
an extra particle on the parent grid and 64 extra particles 
on the sub-grid to serve as seed masses to start the infall. 
The parent seed mass and the sub-grid seed masses all had 
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Figure 4. A phase space plot of the parent grid, and sub-grid 
particles. The open squares are the parent grid particles while the 
open circles are the sub-grid particles. This model was evolved up 
to the onset of non-linearity, S rms = 1.0. 



Figure 6. A comparison of the computed density profiles from 
the nested grid code to the analytic solution. The solid line is the 
analytic solution while the open circles are the sub-grid solution, 
and the filled in boxes are the parent grid solution. 




Figure 5. A plot of the evolved density field on the parent grid 
and the sub-grid at the onset of non-linearity in a ID collapse 
model. No discernible noise can be detected in the sub-grid den- 
sity due to my method of laying down sub-grid particles. The 
filled in circles are the sub-grid densities, while the open boxes 
are the parent grid solution. 



the same coordinates to guarantee that the infall was sym- 
metric with respect to the same point on both grids. After 
several time steps with this perturbed density field the extra 
seed masses were removed and the particles were allowed 
to move based upon the CIC derived density field and the 
resulting potential. The sub-grid was chosen to be 5 par- 
ent grid cells in size with 25 3 sub-grid cells. This implies 
that I can expect about a factor of 5 increase in resolution. 



The amplitude is fixed uniquely by specifying the radius at 
turnaround, Peebles et al. (1989). Furthermore they were 
able to show given an initial mass perturbation mo then the 
turnaround radius, r ta is given by 



m 



— (6tt) p b r ta 



(6) 



Therefore once r ta has been determined, then the dimen- 
sionless radius, A = r/r ta , can be specified and the normal- 
ization of the density profile is set. For this test rta ~ 3.5 
parent grid cells. 

To determine the density profile from the particle po- 
sitions I bin the particles into spherical shells, this gives 
the number density of particles as a function of radius. In 
the case of the parent particles this is the density profile 
(p(r) — p) /p, but in the case of the sub-grid particles this is 
just the number density at radius A. To get (p(r) — p)/p for 
the sub-grid I must divide the number density by the num- 
ber of sub-grid particles per parent particle, this corrects 
for the nl g sub-grid particles which comprise each parent 
particle. This result gives us the final density distribution. 

Fig. 6 is a plot of the evolved density profile for the 
spherical infall problem. In Fig. 7 and 8 I show a dot plot 
of the particle positions within the sub-grid region. Fig. 7 
shows the particle locations on the parent grid for a slice 
2 parent cell in thickness, while figure 7 shows the parti- 
cle positions on the sub-grid 1/4 a parent cell in thickness. 
In each dot plot one can clearly see spherical shells. The 
sub-grid solution has some slight deviations from perfect 
spherical symmetry near the central peak. These deviations 
could be due to my mass advection scheme, and the fact 
that as the particles enter the sub-grid region the interpola- 
tion used to assign sub-grid particle velocities is in error by 
a small amounts. Thus the particles are moving in slightly 
the wrong direction. The effect is small as evidence by the 
accuracy of the density profile, to be considered in Fig 6. 
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Figure 7. A slice plot showing the parent particle positions on 
the sub-grid region for a slice two parent cell in thickness. 



Figure 8. A slice plot showing the sub-grid particle positions 
on the sub-grid region for a slice 1/4 parent cell in thickness. 



The solid line is the Bertschinger analytic solution, the open 
squares are the density profile of the parent particles, and 
the crosses are the sub-grid density profile. There is some 
deviation from the analytic solution for the parent particles 
which is expected since I am using a relatively coarse grid 
for this calculation. The sub-grid appears to agree favorably 
with the analytic solution down to values of A ~ 0.007, which 
corresponds to roughly a length of one-sixth of a sub-grid 
cell. Thus, I appear to be obtaining better resolution than 
one would naively expect from a particle-mesh type code, 
where I would expect the sub-grid solution to be valid only 
down to scales of one sub-grid cell over r ta , or 0.07. In ad- 
dition, the sub-grid solution appears to be roughly a factor 
of 10 better than the parent grid solution. Indicating that 
I am getting better resolution on the sub-grid by roughly a 
factor of 10, when naively based upon the cell sizes on the 
parent grid and sub-grid I would only expect a factor of 4 
in improvement. 



4.4 Spherical Hole 

Another test of how well the code is able to evolve systems 
whose inherent symmetry differs substantially from cubic 
is the spherical hole solution due to Fillmore & Goldreich 
1984), and Bertschinger (1985). For this test I run the code 
with a single sub-grid centered on an uncompensated hole. 
Bertschinger (1985) showed that the central void should 
grow as R oc a™, where n = 1/5. I ran simulations using 
a 16 3 grid with the same number of parent particles on the 
grid. The parent particles and sub- grid particles were then 
binned into spherical shells to compute the density profiles. 
From the profiles I find n = 0.182 ± 0.022, consistent with 
the theoretical result of Bertschinger. 



4.5 Two-Dimensional Perturbations 

For the final test of the code I compare the results of a 512 2 
two-dimensional simulation to the evolution of the same ini- 
tial conditions using my 3D nested-grid code. For this test 
one of the 2D models from Beacom, et al. (1991) was used. 
This model has an initial power-law density fluctuation spec- 
tra, with a spectral index of 0, and a high frequency cutoff 
at k = 32. This choice of cutoff was made so that initially 
there would not be any power in the Fourier modes which 
can only be felt on the sub-grid. Thus, any clustering on the 
sub-grid region will be purely a consequence of evolution. 

The following test is different from the previous test 
cases. In all the previous tests there exists an analytical so- 
lution to the test. Unfortunately all of the analytic solutions 
possess certain symmetries. It is likely that a real model 
will not possess many of the symmetries which have been 
discussed to this point. As a consequence if I wish to com- 
pare solutions to my nested-grid code to more realistic cases 
I must compare to 2D solutions rather than full 3D. This is 
due to the lack of dynamical range in the 3D simulations. 
The 2D code used here for comparisons has been well tested 
in the work of Beacom, et al. (1991), and provides a useful 
benchmark for more sophisticated testing. 

Based upon the final evolved distribution of particles in 
the 2D model, I identify a sub-grid region. Then all of the 
particles which enter the sub-grid are tagged. These parti- 
cles will become the initial distribution of sub-grid particles. 
The full 512 2 distribution is then sampled every fourth par- 
ticle to reduce it to a 64 2 set of particles. These become the 
parent grid particles. The 64 2 distribution of parent particles 
are laid down onto a 64 3 grid. I then define a sub-grid re- 
gion which has the same resolution as the initial 512 2 mesh. 
The nested-grid code is then run up to the same non-linear 
wavelength, in this case k n i = 4, as extrapolated from linear 
theory. 

I can then compare the distribution of particles in the 
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sub-grid region of the original 2D data set to the evolved 
distribution of sub-grid particles on the sub-grid region of 
the full 3D simulation. I should find that the two distribu- 
tions agree nicely. Furthermore, I can use this as a test of 
collisionality, because as particles accrete onto one of the 
clumps in the sub-grid they should not scatter out of the 
plane. Hence in the final 3D distribution none of the sub- 
grid particles should have acquired a velocity component in 
the z direction. 

Fig. 9 is a 64 2 sub-sample of the full 512 2 final evolved 
particle distribution for the original 2D perturbations. The 
left hand column contains the original 2D particles while the 
right hand column contains the particles from the nested- 
grid 3D simulation. The top left is the evolved 2D distri- 
bution plotting every 8th particle. The top right is a slice 
through the parent grid plotting only the sub-grid region. 
The bottom left is the full set of 2D particles on the sub- 
grid region, while the lower right is the evolved sub-grid 
distribution. The bottom two dot pictures of the sub-grid 
region agree nicely with only a slight indication of edge ef- 
fects. Furthermore by comparing the plot in the upper right 
to the plot in the lower right hand corner one can get a 
perspective on the effect of increasing the mass resolution. 

I compare the initial z positions of the particles to the 
evolved z positions of the particles by computing the mean 
of the absolute value of their difference. I find the motion in 
the z direction to of order l/250t/i of a sub-grid cell, and thus 
is completely negligible. There is no motion on the parent 
grid down to round-off level. This indicates that the quasi- 
isolated boundary conditions I apply to the sub-grid region 
introduce little noise in the form of sub-grid particle motion 
out of the plane of the initial 2D perturbations. 

To make the comparison between the dot plots in Fig. 
9 more qualitative I compute the correlation coefficient be- 
tween the density field on the 2d sub-grid region and the 
density field on the 3D sub-grid. To generate a 3D density 
field from the 2D particle positions I extract all particles 
on the sub-grid region plus a small buffer to minimize edge 
effects and stack the appropriate number vertically to make 
the density field three-dimensional. The correlation coeffi- 
cient is defined to be 

K = < 8p2D(x)8psa(x) > 

C2DCTSG 

where 5p2D is the density field generated from the 2D distri- 
bution, Spsa is the density field generated from the sub-grid 
distribution of particles, and oid, and asa are the standard 
deviations of the respective density fields. For the model 
considered here I find K = 0.98. Thus I have quite good 
agreement between the original 2D density field and the full 
3D sub-grid density field. 

4.6 3D Mass Advection Test 

Based upon the ID Zel'dovich test it is clear that the sub- 
grid particles are being advected into the sub-grid region 
correctly when there is little clustering in the sub-grid re- 
gion. Now I wish to test if the same conclusion holds when 
there is clustering of the parent grid particles prior to the 
sub-grid particles being created on the sub-grid. There is to 
my knowledge no analytic solution for such a problem. At 
the request of the referee I am providing this additional test 



of the advection scheme which I am using for the sub-grid 
particles. This test has no known solution and so I will be 
looking for any significant deviations which might be intro- 
duced by the advection scheme. As I will demonstrate I am 
unable to find any such effects. 

To begin the test I use a 32 3 parent mesh with the 
same number of particles. The initial spectrum on the parent 
grid is scale-free and has a spectral index of -1. I evolve the 
system to a non-linear wavenumber of 2 performing data 
dumps at k n i = 16, 8, 4, where k is measured in units of the 
fundamental. Then based upon the k„i = 4 model I find a 
region where a merger has recently taken place between two 
clumps. I define a sub-grid region there and evolve the parent 
grid particles and the sub-grid particles, again perfomring 
data dumps at the same non-linear wavenumbers. In figure 
10 I have plotted a series of dot plots through the sub-grid 
region. The right-hand column is the sub-grid grid particles 
and the left-hand plots are the parent particles. The top 
row are the particle positions at k n i — 16, the middle row 
at k n i = 8, and the bottom row is at k n i = 4. In all cases 
the sub-grid distributions appear to be consistent with the 
parent grid particles. 

Once again to make the arguement more qualitative I 
compute the density cross-correlation function between the 
sub-grid particles and the parent particles in several ways. 
First I extract all the parent and sub-grid particles within 
the sub-grid region plus a small buffer to avoid edge effects. 
I bin both sets of particles onto a mesh the same resolu- 
tion as the parent mesh. I then cross-correlate these den- 
sity fields. These results are in column 2 of table 1, and 
are labeled Ksa- The third column of table 1 contains the 
cross-correlations cimputed in a slightly different manner. 
Since I know from which parent particle each sub-grid par- 
ticle came, I can compute the center of mass of the cloud of 
sub-grid particles for each parent particle on the sub-grid. I 
then bin these centers of mass onto a mesh with the same 
resolution as the parent mesh. I then compute the cross- 
correlation between this new density field obtained from the 
sub-grid particle centers of mass and the parent grid parti- 
cles as before. These cross-correlations are in column 3 of 
table 1, and are labeled Kqm- 

There is relatively good agreement in all cases. Of 
course I don't expect the cross-correlation to be exact. The 
sub-grid has much higher resolution than the parent grid so 
the sub-grid particles will have small scale motions which 
will cause deviations in their motions from what one would 
expect without higher force resolution. Furthermore, the 
agreement between the sub-grid distribution and the parent 
grid distribution gets poorer as the models evolve further. 
This is to be expected due the presence of small scale mo- 
tions on the sub-grid which cannot be resolved on the parent 
grid. In addition, the second set of cross-correlations is bet- 
ter than the first which again is to be expected since the 
center of mass of a cloud of sub-grid particles should remain 
reasonably close to the location of the corresponding par- 
ent particle. Any deviations will be due to the small scale 
motions which can exist on the sub-grid and not the parent 
grid. All of this is consistent with the results of Melott & 
Shandarin (1993), who find that if one computes the cross- 
correlation between models whose only difference is the cut- 
off scale k c then in the limit that k n i « k c the difference 
between the two models should become negligible. 
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Figure 9. A dot plot of a slice through the distribution of particles in the 2D perturbations test. The left hand column contains the 
original 2D distribution, and the right hand column contains the 3D distribution from the nested-grid code. The upper left plot is the 
evolved distribution of particles in the 2D case, where only a 64 2 subset of the original 512 2 particles has been plotted. The small box 
labels the position of the sub-grid region with the parent grid. The upper right is a slice through the sub-grid region of the parent grid 
particles. The lower left is the evolved distribution of particles from the 2D distribution on the sub-grid region. The lower right is the 
evolved distribution of sub-grid particles. The small boxes in the upper left plot labels the sub-grid region. 
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Figure 10. The left-hand column contains plots of the parent particles in a slice 2 parent grid cells thick. The right hand column is the 
sub-grid particle positions in a slice 1 parent cell in thickness. The top row are the particle positions at k nt = 16, the middle row arc the 
particles at fc n ; = 8, and the bottom row are the particles at fc„; = 4. 
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Table 1. Cross-Correlations from the 3D Mass Advcction Test 
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